(* ::Package:: *)

(* ::Section:: *)
(*vertex*)


(* ::Input:: *)
(*ClearAll[*)
(*x1,y1,*)
(*x2,y2,*)
(*x3,y3*)
(*]*)


(*\:5b9a\:4e49\:4e09\:89d2\:5f62\:7684\:4e09\:4e2a\:9876\:70b9*)
x1=-1;y1=0;
x2=2;y2=0;
x3=0;y3=2;


(* ::Section:: *)
(*Triangle element*)


(*\:4e09\:89d2\:5f62\:5355\:5143\:4e2d\:ff0c\:5e38\:7528\:7cfb\:6570\:7684\:5b9a\:4e49---------------*)
a1:=x2 y3-y2 x3;
a2:=x3 y1-y3 x1;
a3:=x1 y2-y1 x2;
(*------------*)
b1:=y2-y3;
b2:=y3-y1;
b3:=y1-y2;
(*------------*)
c1:=x3-x2;
c2:=x1-x3;
c3:=x2-x1;
(*\:4e09\:89d2\:5f62\:7684\:9762\:79ef-----------------*)
\[CapitalDelta]e:=1/2 (b1 c2-b2 c1);


(*\:5750\:6807\:53d8\:6362--------------*)
coToxy[x_,y_]:=1/(2\[CapitalDelta]e)({{b1,c1},{b2,c2}} . {x,y}+{a1,a2});
coTo\[Xi]\[Eta][\[Xi]_,\[Eta]_]:={{c2,-c1},{-b2,b1}} . {\[Xi],\[Eta]}+{x3,y3};


(*test*)
{\[CapitalDelta]e,
coToxy[x,y],
coTo\[Xi]\[Eta][\[Xi],\[Eta]]
}//Expand//Column


Region[
InverseTransformedRegion[
Triangle[{
{x1,y1},{x2,y2},{x3,y3}
}],
(*\:53d8\:6362\:5173\:7cfb*)
(coTo\[Xi]\[Eta][\[Xi],\[Eta]]/.{\[Xi]->Indexed[#,1],\[Eta]->Indexed[#,2]})&
],
Axes->True,
ImageSize->Small
]


Region[
TransformedRegion[
Triangle[{
{x1,y1},{x2,y2},{x3,y3}
}],
(*\:53d8\:6362\:5173\:7cfb*)
(coToxy[\[Xi],\[Eta]]/.{\[Xi]->Indexed[#,1],\[Eta]->Indexed[#,2]})&
],
Axes->True,
ImageSize->Small
]


(* ::Chapter:: *)
(*Gauss quadrature*)


(*\:9ad8\:65af\:79ef\:5206\:516c\:5f0f 1----------------*)
gaussSum[1][fn_]:=Module[{pos,coe},
(*\:91c7\:6837\:70b9-----------*)
pos={
{1/3,1/3}
};
(*\:91c7\:6837\:70b9\:7684\:6743\:91cd------------------*)
coe={1};
(*\:5bf9\:91c7\:6837\:70b9\:5e26\:6743\:6c42\:548c-----------*)
Plus@@(fn@@@pos . coe)
]


(*\:9ad8\:65af\:79ef\:5206\:516c\:5f0f 3 -------------*)
gaussSum[3][fn_]:=Module[{pos,coe},
(*\:91c7\:6837\:70b9-----------*)
pos={
{0.5,0.5},
{0,0.5},
{0.5,0}
};
(*\:91c7\:6837\:70b9\:7684\:6743\:91cd------------------*)
coe={
1/3,
1/3,
1/3};
(*\:5bf9\:91c7\:6837\:70b9\:5e26\:6743\:6c42\:548c-----------*)
Plus@@(fn@@@pos . coe)
]


(*\:9ad8\:65af\:79ef\:5206\:516c\:5f0f 4-----------------*)
gaussSum[4][fn_]:=Module[{pos,coe},
(*\:91c7\:6837\:70b9-----------*)
pos={
{1/3,1/3},
{0.73333333,0.13333333},
{0.13333333,0.73333333},
{0.13333333,0.13333333}
};
(*\:91c7\:6837\:70b9\:7684\:6743\:91cd------------------*)
coe={
-0.56250000,
0.52083333,
0.52083333,
0.52083333
};
(*\:5bf9\:91c7\:6837\:70b9\:5e26\:6743\:6c42\:548c-----------*)
Plus@@(fn@@@pos . coe)
]


gaussSum[7][fn_]:=Module[{pos,coe},
(*\:91c7\:6837\:70b9-----------*)
pos={
{1/3,1/3},
{0.5,0.5},
{0,0.5},
{0.5,0},
{1,0},
{0,1},
{0,0}
};
(*\:91c7\:6837\:70b9\:7684\:6743\:91cd------------------*)
coe={
0.45,
0.13333333,
0.13333333,
0.13333333,
0.05,
0.05,
0.05
};
(*\:5bf9\:91c7\:6837\:70b9\:5e26\:6743\:6c42\:548c-----------*)
Plus@@(fn@@@pos . coe)
]


(*\:9ad8\:65af\:79ef\:5206\:516c\:5f0f--------------------*)
gaussSum[7,2][fn_]:=Module[{pos,coe},
(*\:91c7\:6837\:70b9-----------*)
pos={
{1/3,1/3},
{0.05961587,0.47014206},
{0.47014206,0.05961587},

{0.47014206,0.47014206},
{0.79742699,0.10128651},
{0.10128651,0.79742699},
{0.10128651,0.10128651}
};
(*\:91c7\:6837\:70b9\:7684\:6743\:91cd------------------*)
coe={
0.225,
0.13239415,
0.13239415,
0.13239415,

0.12593918,
0.12593918,
0.12593918
};
(*\:5bf9\:91c7\:6837\:70b9\:5e26\:6743\:6c42\:548c-----------*)
Plus@@(fn@@@pos . coe)
]


(* ::Chapter:: *)
(*Test*)


(*\:6d4b\:8bd5\:79ef\:5206\:7684\:51fd\:6570*)
(*tefn[x_,y_]:=E^(x+y)*)
tefn[x_,y_]:=x^2+y^2+8x*y
(*\:8ba1\:7b97\:79ef\:5206\:7684\:7cbe\:786e\:503c---------------*)
NIntegrate[tefn[x,y],
{x,y}\[Element]Triangle[{
{x1,y1},{x2,y2},{x3,y3}
}]
]//N


tefn2[\[Xi]_,\[Eta]_]:=tefn@@coTo\[Xi]\[Eta][\[Xi],\[Eta]]
jacobi:=2\[CapitalDelta]e


(*\[CapitalDelta]e \:662f\:539f x-y \:7a7a\:95f4\:4e0b\:4e09\:89d2\:5143\:7684\:9762\:79ef*)
gaussSum[7][tefn2]*\[CapitalDelta]e


(* ::Chapter:: *)
(*\:5750\:6807\:53d8\:6362*)


(* ::Input:: *)
(*vecX={1,x,y,z};*)
(*vecL={L0,L1,L2,L3};*)
(*mat={*)
(*{1,1,1,1},*)
(*{x1,x2,x3,x4},*)
(*{y1,y2,y3,y4},*)
(*{z1,z2,z3,z4}*)
(*};*)
(*matInv=Inverse[mat];*)
(*equ=MapThread[Equal,{*)
(*vecX,mat . vecL*)
(*}]*)


(* ::Section:: *)
(*inverse*)


(* ::Input:: *)
(*Det[mat]*matInv//Simplify//TableForm*)


(* ::Section:: *)
(*Jacobi*)


(* ::Input:: *)
(*Lzeta={L1,L2,L3};*)
(*Transpose@D[Rest[mat . vecL/.{L0->1-Total@Lzeta}],{Lzeta}]//TableForm*)


(* ::Input:: *)
(**)


(* ::Input:: *)
(*(y4-y2)(z3-z2)-(y3-y2)(z4-z2)//Simplify*)
